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THE  BIOMECHANICS  OF  SPINAL  AND  HEAD  IMPACT: 
PROBLEMS  OF  MATHEMATICAL  SIMULATION 


Y.  King  Liu,  Ph.D. 

Assoc.  Prof,  of  Biomechanics 
Tulane  University  School  of  Medicine 
New  Orleans,  Louisiana  70112 


ABSTRACT 


ii/The  present  expository' paper  examines  the  various  mathematical 
models  proposed  in  impact  studies  in  general  and  those  in  connection 
with  spinal  and  head  injuries  in  particular.  First,  the  concept  of 
injury  tolerance  surface  is  introduced.  When  a  complete  mapping  of 
this  surface  is  in  hand,  then  the  probability  of  injury  due  to  any 
acceleration  vector  can  be  evaluated.  After  reviewing  a  typical  single- 
degree-of-freedom  model,  a  quadrature  scheme,  based  on  classical  Fourier 
transform  technique,  is  proposed  for  obtaining  the  pulse  response  from 
the  experimentally  determined  mechanical  impedance.  The  ignoring  of 
the  spatial  mass  distribution  is  shown  to  be  the  source  of  many  dif¬ 
ficulties  in  simple  models.  Typical  results  for  the  pilot-ejection 
problem  where  the  spatial  mass  distribution  is  accounted  for  is  then 
given.  The  determination  of  input  parameters  necessary  to  implement 
this  class  of  models  is  also  briefly  described. 

The  second  part  of  the  paper  begins  with  a  review  and  critique 
of  a  promising  one-dimensional  continuum  model  of  head  injury  by  Hayashi*4. 
Exact  wave-propagation  solutions  were  obtained  for  the  intracranial  pressure 
and  container  acceleration.  When  compared  with  the  approximate  solution 
given  previously,  both  mathematical  and  physical  reasons  were  advanced  for 
the  improvement  of  this  model.  The  relatively  small  number  of  dimensionless 
parameters  of  the  anticipated  improved  Hayashi  model  will  make  any  future 
discussions  of  the  cavitation  theory  very  much  easier.  The  survey  of 
two-  or  three-dimensional  models  is  divided  into  axisymmetric,  rotational 
and  nonaxisymraetric  versions  of  fluid-filled  spherical  shells.  A  brief 
summary  of  the  implications  of  the  idealizations  discussed  above  to  other 
critical  organs  of  the  human  body  concludes  the  paper. 


INTRODUCTION 


In  the  course  of  daily  living,  a  host  of  circumstances  confront 
the  human  body  such  that  a  mechanical  input  is  imposed  on  it.  Many 
are  accidental,  e.g.,  athletic  injuries,  automobile  and  industrial 
accidents,  while  others  are  preplanned,  e.g.,  ejection  from  high-speed 
aircraft,  manned  launchings  and  recoveries  of  aircraft  and  spaceships 
on  water  or  land.  In  either  case,  the  knowledge  of  the  response,  toler¬ 
ance  and  control  of  the  body  is  highly  desirable,  and  in  many  situations, 
quite  imperative.  In  terms  of  its  response  to  externally  applied  forces, 
the  body  obeys  the  laws  of  mechanics.  There  is  little  doubt  that  this 
particular  structure  will  be  with  us  for  a  long  time  to  come.  Any  un¬ 
derstanding  or  knowledge  gained  with  respect  to  its  structural  mechanics 
will  have  a  continuing  significance.  The  analogy  that  the  skeletal 
(bone)  system  corresponds  to  the  main  structural  members  of  machines  is 
a  plausible  one.  In  that  spirit,  ligaments  and  tendons  are  then  cables 
and  tie-bars  and  muscles  are  sources  of  motive  power.  It  is  inappropri¬ 
ate  to  carry  the  analogy  too  far.  Examples  of  physiological  function 
which  has  rare  counterparts  in  machines  are:  (1)  The  neuromusculature 
can  shorten  as  much  as  60%  of  its  resting  length  due  to  active  contrac¬ 
tion  and  generally  sustain  loads  several  times  more  than  its  bony  part. 
(2)  Living  bone  can  repair  itself  after  fracture. 

The  basic  modeling  objectives  are  two-fold:  (a)  To  propose  the 
simplest  mathematical  description  consistent  with  explaining  the  physi¬ 
cal  phenomenon.  At  the  core  of  the  mathematical  description,  there  is 
the  anatomical  and/or  physiological  framework  upon  which  the  model  rests 
and,  on  the  other  hand,  the  biomechanical  data  needed  to  implement  the 
model,  (b)  To  provide,  through  the  model  or  its  improvement,  new  in¬ 
sights  or  predict  new  phenomena,  which  serve  to  indicate  where  the 
measuring  devices  should  be  placed.  The  development  of  meaningful  ex¬ 
periments  should  be  guided  by  some  theory,  no  matter  how  simple.  The 
construction  of  a  meaningful  theory,  on  the  other  hand,  requires  the 
availability  of  experimental  data  and  clinical  observations.  If  data 
are  not  provided  or  found,  there  is  always  the  danger  that  the  theory 
will  be  based  on  how  the  body  could  or  should  respond  rather  than  on 
how  it  does. 


INJURY  TOLERANCE  SURFACES 

The  mathematical  models  of  the  human  body  subjected  to  accelera¬ 
tions  depend  on  the  direction  of  the  acceleration  vector.  The  applied 
accelerations  range  from  steady-state  oscillations  to  abrupt  pulses.  For 
the  sake  of  uniformity  in  the  discussion  to  follow,  we  shall  follow  the 
physiological  acceleration  coordinate  (G)  system  recommended  by  the 
Committee  on  Acceleration  of  the  AGARD  Aerospace  Medical  Panel  as  report¬ 
ed  by  Gell*. 

The  fundamental  problem  in  acceleration  injury  is  the  construction 
of  tolerance  surfaces  in  the  3-G  space  with  a  suitable  quantity  as  pa¬ 
rameter.  For  abrupt  pulses,  the  time  of  uniform  acceleration  exposure 
is  an  obvious  candidate  as  a  parameter.  This  point  of  view  is  easily 

illustrated  in  the  G  -G  plane.  Coincidentally,  most  of  the  current 
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available  laboratory  data  on  acceleration  injury  tolerance  are  in  these 
directions  of  impact.  In  Figure  1  the  data  of  Eiband^  is  replotted  for 
a  uniform  acceleration  duration  of  0.05  secs.  The  axes  intercepts  cor¬ 
respond  to  the  known  tolerance  data  for  ±GZ.  The  straight  dotted  lines 
represent  the  easiest  analytical  construction  of  this  tolerance  curve 
and  is  in  no  way  suggestive  of  reality.  In  fact,  from  the  F4C  aircraft 
data  in  which  the  ejection  seat  is  inclined  14°  towards  the  +x  axis,  a 
second-degree  hyperbola,  shown  as  a  solid  curve  might  be  much  more  ap¬ 
propriate.  Tolerance  data  for  combinations  of  accelerations  are  woe¬ 
fully  lacking  even  for  the  plane  and  is  nonexistent  for  the  3-G  space. 
Many  acceleration  environments  are  indeed  in  this  category,  e.g.,  (1) 
winged-aircraft  and  automobile  collisions  are  points  in  the  Gx-Gy  plane; 
(2)  helicoptor  crashes  are  points  in  the  Gz-Gx-Gy  space.  This  concept 
of  the  tolerance  surface  is  adapted  from  the  theory  of  yield  surfaces 
(in  the  principal  stress  space)  in  elasticity  and  plasticity. 

REVIEW  OF  SIMPLE  LUMPED- PARAMETER  MODELS 

Given  a  typical  impact  situation,  one  can  usually  discern  the 
following  elements:  mass,  elasticity,  dissipation  and  the  nature  of 
the  applied  impulse,  which  contribute  to  the  dynamic  response  of  the 
body.  Assuming  linearity  of  the  elements,  one  can  obtain  the  response 
of  a  damped  dynamic  system  to  a  given  acceleration  (or  force)  input, 
a(t),  as  the  solution  of  the  well-known  differential  equation: 

x  +  26u0x  +  *  a(t),  (1) 

where  u>0  is  the  natural  frequency,  6  is  the  damping  ratio  and  x(t)  is 
the  xelative  displacement  of  the  support  with  respect  to  the  mass.  The 
solution  of  equation  (1)  for  an  arbitrary  a(t)  is  in  the  form  of  Duha- 
mel's  integral: 

x  =  D  exp  {(-6  +  i  \j  1  -  d2)a)0t}  +  /a(x)h(t  -  x)dx,  (2) 

where  D  is  an  arbitrary  constant  dependent  on  the  initial  conditions, 
i  -  (-1)*5,  exp  =  e  and  h(t)  is  the  indicial  response,  i.e.,  the  response 
to  a  Dirac-delta  input,  of  the  system.  For  quiescent  initial  conditions, 
it  has  been  shown  by  many  investigators,  e.g.,  von  Gierke-*,  Kornhauser4, 
PayneS,  that  the  system  response  can  be  divided  into  approximately  two 
regions  at  a  critical  pulse  duration,  tc.  For  t  <  tc,  influence  on  the 
system  occurs  only  in  the  pulse  area,  which  is  equal  to  the  imposed  vel¬ 
ocity  change.  For  t  >  tc,  the  pulse  shape  is  the  important  factor.  If 
one  were  to  assume  that  the  elastic  element  always  breaks  at  a  given 
peak  force  level,  the  so-called  equal  tissue  strain  assumption,  it  is 
then  possible  to  relate  the  damage  or  injury,  through  the  system  response, 
to  the  parameters  of  the  input  pulse.  Typically,  the  result  is  a  plot 
of  the  maximum  applied  acceleration,  Am,  versus  the  duration  of  the  pulse, 
ti,  on  a  log-log  scale,  as  shown  in  Figure  2.  In  this  tolerance  graph, 
the  "knee"  represents  the  critical  duration,  tc.  To  its  left,  the  re¬ 
gion  is  sensitive  to  velocity-change  and  to  its  right,  the  peak  accel¬ 
eration  of  the  pulse. 

An  alternative  to  the  above  is  the  use  of  response  to  a  sinusoidal 
input  to  predict  system  response  under  arbitrary  pulse  loading.  Let 
the  input  be  the  real  part  of 
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a(t)  =  Ae  (3) 

where  A  is  a  complex  constant.  The  solution  to  (1)  subjected  to  the  in¬ 
put  of  (3)  may  be  written  as: 

x  *  B  exp  {(-6  +  i  ^1  -  S2u>0)t}  +  Aela)t/Z(u) ,  (4) 

where  B  is  a  complex  constant  dependent  on  the  initial  conditions,  Z(u) 
is  the  characteristic  impedance  and 

Z(u)  =  m(u2  -  u2  +  2i6uQu) .  (5) 

The  reciprocal  of  (5),  i.e.,  H(o>)  =  1/Z(a>),  is  the  transfer  function  or 
frequency  response  function.  For  steady-state  motion,  the  transfer  func¬ 
tion  is  the  ratio  of  the  output  (displacement)  to  the  input  (sinusoidal 
excitation) . 

If  one  attempts  to  experimentally  identify  the  transfer  function, 
an  immediate  difficulty  arises.  Because  the  human  body  parts  consist 
of  distributed  masses,  generally  no  point  can  be  picked  where  the  dis- 
placements  ot  the  effective  mass  can  be  measured.  Coermann^  circum¬ 
vented  this  difficulty  by  adopting  the  concept  of  mechanical  impedance. 
Specifically,  mechanical  or  driving  point  impedance  is  defined  as  the 
ratio  of  the  transmitted  force,  Ftr,  to  the  velocity,  x,  of  that  point 
where  the  force  is  transmitted.  It  is  equivalent  to  electrical  imped¬ 
ance  if  one  were  to  choose  the  force-current  analogy.  This  concept  of 
mechanical  impedance  is  not  identical  to  characteristic  impedance,  Z(u), 
defined  in  (4).  To  avoid  confusion,  we  denote  mechanical  impedance  by 
Y(w),  which  is  easily  shown  to  have  the  following  form: 

(1  +  2i6u )/u>  )mu> 

Y(u)  =  - 2 -  (6) 

{1  -  (u/w  )2}  +  2i5u/w 
o  o 


or 


|Y(u)| 


mu 


1  +  (25u/uo)2 

[l  -  (u/uo)2]2  +  (25u/uo)2 


(7) 


The  value  of  frequency  ratio,  u/u0  =  p,  at  which  the  above  equation 
attains  its  maximum  is  found  by  setting  dY/dp  =  0.  It  can  be  shown  to 
occur  at 

p  =  { 46 2  +  (1  +  852)Js}/{1  +  862  -  1664}  .  (8) 

Using  the  above  model,  Coermann^  vibrated  human  volunteers,  dummies  and 
animals  in  the  sitting  and  standing  position  from  1  to  20  hertz.  From 
the  experimentally  obtained  frequency  ratio  corresponding  to  the  peak  of 
the  mechanical  impedance  curve,  he  determined  the  damping  ratio,  6,  with 
the  help  of  (8). 


An  interesting  question  suggests  itself  at  this  point:  given  the 
frequency  response  or  the  mechanical  impedance  of  a  system,  what  can  one 
say  about  its  response  to  an  arbitrary  pulse?  Define 
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a (u)  =  (l/2ir)  I  a(t)e  dt 

J  0 

as  the  Fourier  Transform  of  a(t).  Its  inverse  transform  is: 


■C 


— .  +io)t, 
a(o>)e  dui  . 


Implicit  in  the  above  definitions  is  that  the  conditions  for  the  exis¬ 
tence  of  such  a  transform  pair  are  satisfied  and  that  a(t)  =  0  for 
t  <  0.  The  above  synthesis  when  superposed  on  (4)  yields 


x  *  B  exp  {(-6  *  i  \  1  -  62)u)Qt}  + 


H(uj)a(a))elaitda)  . 


In  words,  the  pulse  response  of  system  (1)  is  the  sum  of  the  transient 
solution  and  the  inverse  Fourier  transform  of  the  product  of  the  fre¬ 
quency  response  function,  H(w)  and  the  frequency  spectrum,  a(io),  of  the 
pulse.  To  compute  the  displacement  response  of  the  effective  mass  to 
an  arbitrary  pulse  based  on  frequency  response  data,  one  needs  only  to 
compare  (5)  and  (6)  and  note  that 

Y(u>)  =  m2ua)2{l  +  2i6  (u>/u>0)  >/Z(uj)  .  (12 

Neglecting  the  transients,  we  get 


■s: 


+“  Y(u)a(u3elwtdt> 


00  m2o)2o){l  +  2i6(oj/w0)} 


For  specific  analytical  forms  of  Y(u)  and  a(w),  it  may  be  possible  to 
obtain  an  analytical  expression  for  x.  In  other  instances,  the  integrand 
in  (13)  is  so  complicated  that  the  evaluation  of  the  integral  is  a  mat¬ 
ter  of  great  difficulty.  The  true  value  of  (13)  in  the  present  problem 
is  to  use  the  experimentally  determined  mechanical  impedance,  Y(w)  to 
yield  information  about  its  response  to  an  arbitrary  pulse.  For  these 
cases,  recourse  generally  has  to  be  made  to  numerical  methods.  To  ac¬ 
complish  the  complex  integration  in  (13),  a  feasible  approach  is  to  re¬ 
solve  the  complex  function  into  separate  functions  of  real  and  imaginary 
part  and  integrate  each  separately.  Rewriting  (13)  as 

i+oo 

F(u))e10)tdui,  (13') 

where  F (to)  =  A(oj)  +  iB(w).  It  follows  therefore,  that 

x (t)  =  a(t)  +  ib(t)  «  (a2(t)  +  b2(t)}^*ei*(t)  ,  (14) 

where 

i+co 

{A(to)  coshwt  -  B(<j)  sinhutjdo) 
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b(t)  - 


<J>  *  tan_1(b/a)  . 

The  practical  difficulty  lies  in  the  frequency  range  to  be  swept.  If 
the  test  data  obtained  is,  say,  from  0  -  200  hertz,  then  the  Hobson's 
choice  is  to  sweep  from  -200  to  +200  hertz  in  evaluating  (14) .  The 
numerical  quadrature  required  can  be  achieved  through  any  number  of 
standard  integration  routines.  The  use  of  a  hybrid  analog-digital  com¬ 
puter  should  also  accomplish  the  task  in  one  operation,  i.e.,  from  the 
data  obtained  to  the  final  tolerance  graph. 

The  massive  experimental  program  mounted  on  behalf  of  the  single- 
degree-of-freedom  model  served  to  identify  some  of  the  critical  organ 
systems  associated  with  a  given  force  environment.  The  best  understood 
case  is  the  caudocephalad  or  +GZ  acceleration  where  the  most  critical 
organ  is  the  spine  and  the  mode  of  injury  is,  with  few  exceptions,  com¬ 
pression  fracture  of  the  vertebral  end-plates.  Given  this  critical  in¬ 
jury  information,  it  is  possible  to  refine  the  model  analysis  and/or 
identification  of  the  model  parameters.  Liu?  recently  reviewed  the 
modelling  and  experimental  progress  made  for  the  spine  in  the  dynamic 
environment  with  particular  emphasis  in  the  +GZ  orientation.  A  brief 
summary  of  his  conclusions  appear  in  the  next  section  of  this  paper. 

The  critical  organs  and  failure  criterion  for  other  orientations  of 
impact  are  only  partially  defined.  It  is  in  these  orientations  where 
the  relationships  given  by  equations  (2)  and  (13)  are  still  needed  pend¬ 
ing  additional  work  which  will  delineate  the  failure  modes  in  the  criti¬ 
cal  organs.  As  an  example,  consider  the  dynamics  of  the  head-neck  junc¬ 
tion,  which  has  major  implications  in  brain  injury  due  either  to  whip¬ 
lash  or  a  direct  blow  to  the  head.  Initially,  it  can  be  viewed  as  a 
single  degree  of  freedom  rotational  problem,  i.e.,  the  dependent  vari¬ 
able  in  (1)  is  the  angular  displacement,  0,  of  the  head.  The  driving 
point  impedance  has  been  determined  in  detail  by  Hodgson  et  al.8  As¬ 
suming  a  linear,  single-degree-of-freedom  system  the  parameters,  6  and 
uj0,  can  be  identified  from  the  impedance  data  and  then  the  response  to 
a  pulse  is  found  from  (2).  Alternately,  if  one  were  to  compare  the  ex¬ 
perimentally  determined  impedance  data  with  the  response  as  given  in  (4), 
the  data  will  indicate  its  deviations  from  the  simple  model.  These  devi¬ 
ations  can  either  be  due  to  additional  degrees  of  freedom  in  the  system 
and/or  nonlinear  material  properties  or  geometry. 

MODELS  CONSIDERING  MASS  DISTRIBUTION 

Liu?  recently  gave  a  review  and  an  assessment  of  the  multi-degree- 
of-freedom,  discrete-parameter  and  continuous-parameter  models  of  the 
spine  under  inertial  loading  and  their  relationships  to  the  simpler 
models  proposed  earlier.  The  essential  points  made  in  that  paper  were: 

(1)  Most  of  the  inconsistencies  of  the  simple  models  were  due  to 
ignoring  the  spatial  mass  distribution. 

(2)  The  possible  occurance  of  shock  waves  in  a  lightly  damped  non¬ 
linear  rod  model  of  the  spine  is  due  to  an  accumulation  of  small  effects. 


{A (m)  sinhut  +  B(<d)  coshu>t}du> 
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The  shock  forms  as  a  result  of  the  steepening  of  the  negatively  sloped 
portion  of  the  input  pulse  and  not  as  a  function  of  rapid  rise  time 
alone  as  has  been  suggested  in  the  literature. 

(3)  That  a  lumped-parameter  equivalent  to  a  continuous  segment  is 
valid  provided  its  length  is  much  shorter  than  the  shortest  wavelength 
of  interest.  If  one  were  to  have  an  understanding  of  the  spinal  injury 
mechanism,  the  wavelength  of  interest  is,  at  least,  the  thickness  of  the 
disc. 

(4)  The  initial  curvature  and  the  eccentric  center  of  mass  of  the 
trunk  (anterior  to  the  vertebral  bodies)  places  the  spine  under  complex 
loading  even  under  normal  conditions.  For  impact  loading,  the  situation 
is  greatly  exaggerated. 

(5)  In  order  to  implement  either  the  multi-degree-of-freedom  con¬ 
figuration  model  of  Ome  and  Liu9  or  the  distributed-parameter  model  of 
Moffatt*0,  a  large  amount  of  biomechanical  data  is  needed,  e.g.,  the 
inertial -property  distribution  of  the  head,  neck  and  trunk,  the  failure 
surface  of  the  intervertebral  joint  under  complex  loading  and  the  initial 
configurations  of  the  pilot  prior  to  ejection,  etc.  Liu  et  al.H  recent¬ 
ly  obtained  the  inertial  properties  of  a  segmented  cadaver  trunk.  Fig¬ 
ures  3-5  summarized  their  results  for  the  distribution  of  mass,  the 
centers  of  mass  and  moments  of  inertia  respectively. 

(6)  The  predominance  of  anterior- lip  and/or  compression  fractures 
in  the  lower  thoracic  vertebrae  during  pilot  ejection  are  primarily  at¬ 
tributable  to  the  large  negative  moments  there  and  only  secondarily  to 
the  axial  and  shear  forces.  These  summary  results  are  shown  in  Figures 
6-8.  The  solid  curves  are  the  results  obtained  by  using  the  inertial 
data  of  Liu  et  ul.H  as  input  to  the  model  by  Ome  and  Liu9  and  the 
dotted  curves  were  the  corresponding  results  using  assumed  data.  Figure 
9  illustrates  the  time  history  of  the  configuration  changes  of  the 
spine  under  these  same  loading  conditions. 

HEAD  IMPACT  MODELS 

The  special  vulnerability  of  the  head  to  injury  as  compared  to  other 
parts  of  the  human  body  is  evidenced  by  the  fact  that  about  75%  of  the 
fatalities  from  all  accidents  are  linked  to  craniocerebral  trauma.  The 
Proceedings  of  the  Conference  on  Head  Injury,  edited  by  Caveness  and 
Walker^  summarized  the  spectrum  of  activities  in  connection  with  the 
problem  up  to  1966.  Recently,  GoldsmithlS  reviewed  the  more  current 
work.  It  would  be  superfluous  to  cite  these  contributions  once  more 
except  perhaps  in  illustrating  some  inconsistencies  which  have  entered 
the  modeling  of  intracranial  pressure  in  head  impact.  The  main  discus¬ 
sion  will  deal  with  models  which  have  entered  the  literature  since  1966. 

Mathematical  models  of  the  head  have  been  divided  in  much  the  same 
fashion  as  the  major  hypotheses  proposed  to  explain  head  injuries  due 
to  impact.  One  group,  the  rotational  school,  has  contended  that  the 
rotational  acceleration  induced  by  impact  causes  high  shear  strains  in 
the  brain  matter,  rupturing  cerebral  blood  vessels  and  tissue.  The 
cavitation  school,  on  the  other  hand,  has  claimed  that  there  exist  points 
within  the  brain  where  the  reduced  pressure  is  sufficient  to  rupture 
the  capillary  walls.  The  normal  transmural  pressure  of  these  capillaries 
is  a  few  mm  Hg;  however,  due  to  the  impact,  the  transmural  pressure  is 
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suddenly  Increased,  bursting  the  capillary.  Of  course,  when  the  pressure 
is  reduced  to  near  vapor  pressure  of  the  brain  substance,  cavitation 
should  take  place.  The  catastrophic  collapse  of  the  bubbles  thus  formed 
is  another  possible  cause  of  brain  damage. 

ONE  DIMENSIONAL  TAVITATION  MODEL 


The  most  recent  one-dimensional  continuum  model  was  given  by  T. 
Hayashi^^,  The  system  consists  of  a  rigid  but  massless  vessel  (skull) 
containing  elastic  fluid  (brain  and  cerebrospinal  fluid) .  The  vessel  is 
attached  to  a  linear  spring  k,  which  represents  the  composite  elastic 
properties  of  the  helmet,  skull,  hair  and  elasticity  of  the  wall.  Thus, 
the  problem  can  be  simplified  to  that  of  a  fluid  "rod"  enclosed  in  a 
rigid  vessel  with  an  attached  spring  striking  a  rigid  wall,  see  Figure 
10.  The  governing  nondimensional  differential  equation  is: 


u  +  £ 

TT  TT 


’XX* 


(16) 


where  the  subscript  is  used  to  denote  partial  differentiation  with  re¬ 
spect  to  the  independent  variables,  x  and  t.  The  associated  initial  and 
boundary  conditions  are  respectively: 

C(x,0)  =  u(0)  =  0;  ?T(x,0)  -  0;  ut(0)  =  V 

C(0,t)  -  C(1,t)  =  0;  u(t)  =  AUx(0,t)  -  Cx(1,t)}  . 

In  terms  of  the  original  quantities  (with  tilde  ~) ,  the  variables 
(18)  have  been  nondimens ionali zed  as: 

?  =  £/A;  u  =  u/A;  x  =  x/A;  t  =  cT/A 

c  =  (B/p)5*;  kf  =  BA/A;  A  »  kf/k;  V  =  vQ/c, 

where  ^(xYr)  is  the  displacement  of  the  fluid  at  location  K  relative  to 
the  vessel,  A  is  the  area,  A  is  the  length,  T  is  the  time,  TI  is  the  abso¬ 
lute  rigid-body  displacement  of  the  vessel,  v0  *  velocity  of  the  vessel 
just  prior  to  impact;  B  is  the  bulk  modulus;  p  is  the  density;  c  =  (B/p)5* 
is  the  wave  speed  and  k£  is  the  stiffness  of  the  fluid.  Hayashi*4  utiliz¬ 
ed  a  separation  of  variables  technique  to  obtain  an  infinite  series  solu¬ 
tion  for  the  vessel  displacement  and  the  fluid  pressure  in  the  original 
variables,  i.e.,  u  and  In  terms  of  the  dimensionless  variables  given 

in  (19),  his  solutions  can  be  written  as: 

,»  sin(2u)  ) 

u(t)  =  2V  V  - - - sin(2o>nT)  (20) 

n*l  (2u)n){2u)n  +  sin(2un)} 


(17) 

(18) 

in  (16)- 
(19) 


ana 


-P(x,t) 


CX(X,T) 


■A  sin(o  ) 

4V  S  2„  ♦  sinft»  )  (l-2x)}  sin(2V).  (21) 

n=l  n  v 
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where  w  are  the  n  roots  of  the 
n 


characteristic  equation: 


a)  tamo 


_1 

4X 


(22) 


The  roots  of  (22)  are  partially  tabulated  in  the  National  Bureau  of  Standards 
Handbyok  of  Mathematical  Functions  (edited  by  Abramowitz  and  Stegun).  Any 
attempt  to  perform  computations  using  (20)  -  (22),  however,  would  show 
that  these  infinite  series  solutions  are  slowly  convergent,  if  at  all, 
especially  for  small  values  of  time.  In  short,  while  (20)  and  (21)  are 
nominally  "exact"  solutions,  they  are  useless  in  computing  numerical 
results.  If,  in  addition,  certain  approximations  are  made  in  order  to 
make  (20)  -  (22)  more  pliable,  the  errors  are  then  compounded  as  will  be 
shown  presently  by  the  exact  closed  form  solution  obtained  by  the  author. 


Taking  the  Laplace  transform  of  Equations  (16)  to  (18)  yields 


r(x,p) 


V{cosh  px  -  tanh(p/2)  sinh  px  -  1} 
p2 { 1  +  2Ap  tanh(p/2)} 


and 


(23) 


u(p) 


2XV  tanh(p/2) 
p(l  +  2Ap  tanh(p/2)} 


(24) 


The  pressure  in  the  fluid  is  found  from  the  inversion  of: 

*«.p>  -  ■  vlpu .%  s$fjrh  p’i)  • 


(25) 


The  poles  of  (25)  are  the  characteristic  roots  of  (22) .  Using  the  standard 
contour  inversion  integral  and  summing  the  residues  will  yield  identical 
results  as  those  shown  in  (20)  and  (21) .  We  are  back  in  the  same  quandary! 
Under  these  circumstances,  instead  of  dealing  with  the  entire  function,  one 
circumvents  the  problem  by  dealing  with  its  expansion,  see  Carslaw  and 
JaegerlS.  Expanding  the  hyperbolic  functions  in  terms  of  exponentials, 
we  write  (25)  as: 


P(x,p) 


p(l+2Ap){l+e (p)e'P) 


where 


e(p) 


-  1 


1  +  2Ap 

Using  the  binomial  theorem,  we  get 


{l  +  e(p)e'p  i  _1  =  1  + 

'  '  n*l 


(-l)nen(p)e'np 


Observe  that 

c-i)V(p)  •  {i  - }  ■  i  ♦  £  (-1)’ 


0 


(26) 


(27) 


(28) 


(29) 
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where 


O'sn&jr 


are  the  binomial  coefficients. 


In  view  of  (28)  and  (29) ,  (25)  becomes 


-V 


P(x,p)  *  p(l+2Xp) 


it 


e-px_e-p(l-x)- 


fl  ♦  £  e_p(n+x)  -  £  e‘p(n+1 

n*l  n*l 


-x) 


t  i  <p*2r>'“[e'p<n**)-'"p(n'1"’0]  I 


(30) 


—  „  ,  ,  ,  l.v+1  . 

The  inverse  of  l/pCp+jp  is 


C3ePT/p(p^,v*1>p  ■  (2i)V,1h‘t/2>  4]  <3i) 


j  J-Y+l" 
2iri 


Taking  advantage  of  the  shifting  theorem,  the  exact  solution  to  the  problem 
posed  by  (16)  -  (18)  is: 

.  ZilLLL  *  jl-e"(T“x)/2A]H(T-x)  -  |l-e‘fT‘1+x)/2XjH(x-l+x) 


Si 


l_e-(T-n-x)/2X 


]H(r-n-x)  -^[l-e 


-(x-n-l+x)/2X 


1+x) 


i  V*  Vf  nvA  1  ITi  r-(T-n-x)/2X  y*  (x-n-x)w1  , 

*  v?i(  )  Ulwil  h insrj^ > 

-  re-(T-n-l«)/2X  £  (x-n-l+x)*1 1 H(T-n-l*x) 

L  p-0  p!(2X)V  J 


(32) 


Where  H(t)  is  the  Heaviside  function,  i.e.,  it  is  zero  when  t  <  0  and  unity 
when  t  >  0. 


We  note  that  (32)  is  an  odd  function  about  x  ■  h,  i.e., 

P(x,x)  =  -P(l-x,x)  (33) 

and  in  particular 

P(0,x)  =  -P(l,x).  (33*) 

For  x  =  h,  Pfajx)  =  -P^.t)*  which  is  possible  if  and  only  if  P(%,x)  s  0. 
Equation  (33)  stares  that  the  pressure  at  the  center  is  always  zero  and  the 
pressure  variations  in  time  are  equal,  except  for  sign,  on  either  side  of 
the  center.  Further,  the  solution  (32)  is  valid  only  prior  to  rebound  which 
occurs  when  F  =  kfl  >  0.  In  view  of  (18)  and  (33'),  the  rebound  condition  is 
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equivalently 

-P(0,t)  i  0,  (34) 

i.e.,  when  the  system  returns  to  the  initial  position  of  the  container 
prior  to  impact.  Due  to  the  nature  of  the  Heaviside  function,  for  a  given 
x  and  t,  only  a  finite  number  of  terms  are  needed  from  (32)  to  perform  the 
numerical  calculations,  attesting  to  its  simplicity  as  well  as  exactness. 
For  up  to  2  wave-transit  times  even  hand  calculations  would  suffice. 


The  displacement  can  be  found  immediately  from  the  boundary  condition 
(18),  i.e., 


u(t)  -  Xkx(0,T)  -  ?x(1,t) }  -  -2AP(0,t),  (35) 


where  use  was  made  of  F(x,r)  »  -Cx(x,t)  and  P(0,t)  ■  -P(1,t).  The  differentiation 
in  the  t  domain  corresponds  to  multiplication  by  p  since  the  initial  dis¬ 
placement  is  zero,  hence  ( 

pu(p)  =  2AV  tanh(p/2)(l  +  2Ap  tanh(p/2)}_1  .  (36) 

To  check  that  the  initial  velocity  is  satisfied,  we  note 

u  (0)  ■  lim  p2u(p)  *  lim(2AVp  tanh(p/2) [l+2Ap  tanh(p/2)]  *} 

p-H»  p-H» 


*  lim 

p-H» 


2AV  tanh(p/2) 

1/p  ♦  tanh (p/ 2) 


V. 


(37) 


The  acceleration  is,  therefore,  the  inverse  transform  of 


I(P)  ■  p!u(p)  -  uT(0,  -  (W 

Through  an  expansion  and  inversion  procedure  similar  to  (26)-(29),  we  obtain 

~2XyW  .  e"T/2XH(x)  *  2^)  e_(T'n)/2XH(t-n) 

n*=l 

+  E  E  ~7~  |(T-n)e'(T'n)/2XH(T-n) 

n«l  v®l  \v/  A  v!  1 

-  (T-n-l)e_(T"n'1)/2XH(T-n-l)J  .  (39) 


Equations  (39)  and  (32)  constitute  the  exact  principal  results  of  the 
model.  The  corresponding  Hayashi  approximate  results  are  found  from  (20), 
(21)  and  (22)  by  the  substitution  of  u>j  a  1/2A**  *nto  the  infinite  series 
truncated  after  one  term,  i.e.,  f 
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2Xa 

V 


2\*  sinCi/xS 


(20') 


.  Pfe,T.).  „  lilgO  sin(T/x>i)  (21 ') 

v  2X^ 

Figure  11  is  a  plot  comparing  the  acceleration  of  the  container  as 
obtained  by  equation  (39)  and  (20')  for  different  values  of  stiffness  ratio  X. 
Whereas  the  exact  wave  propagation  solution  has  a  discontinuous  jump  at 
t  =  0  and  "jumps"  at  t  *  1,2,.  .  .,  the  points  of  wave  reflection,  in  contrast, 
the  Hayashi  approximation  is  a  sine  wave.  The  discontinuous  jumps  shown 
in  Figure  11  have  major  implications  for  the  infinite  series  solutions  such 
as  (20)  and  (21).  Would  the  addition  of  more  terms  from  (21)  improve  the 
acceleration  result?  It  is  well-known,  however,  that  the  partial  sums  Sm(x) 
of  a  Fourier  series  for  a  periodic  function  f  cannot  approach  f(x)  uniformly 
over  an  interval  that  contains  a  point  where  f  is  discontinuous.  The  nature 
of  the  deviation  of  Sa(x)  from  f(x)  on  such  intervals  is  known  as  Gibb's 
phenomenon,  see  Carslawl7.  The  pathological  nature  of  the  exact  acceleration 
solution  was  due  to  the  absence  of  the  container  mass,  which  is  to  be 
regarded  as  a  weak  solution  in  the  sense  of  Courant  and  Hilbert1®. 

For  the  pressure  results,  two  sets  of  graphical  data  are  of  interest: 

(a)  the  pressure  field  variation  at  given  times,  and  (b)  the  pressure-time 
variation  at  a  given  location,  x.  In  Figure  12  is  shown  the  dimensionless 
contrecoup  (x  *  0)  pressure  as  a  function  of  t.  When  the  contrecoup  pressure 
reaches  zero,  the  solution  (32)  is  no  longer  valid  because  of  the  rebound 
condition  (34).  Note  that  the  softer  the  impact,  (increasing  X)  the  longer 
the  contact  duration,  tc,  e.g.,  X  *  1,  tc  -  3.275,  whereas  X  *  0.1,  tc  =  1.34. 
For  very  soft  impacts,  X  >>  1,  the  pressure  wave  would  have  traversed  the 

fluid  length  many  times  during  the  time  the  vessel  takes  to  return  to  its 

original  position,  i.e.,  just  prior  to  loss  of  contact.  Also  shown  in  Figure 
12  is  a  comparison  of  the  exact  wave-propagation  solution  with  the  Hayashi 
approximate  solution  for  various  X.  We  note  with  interest  that  for  X  *  10 
the  two  results  are  indistinguishable  and  that  for  X  =  1,  the  approximation 
is  still  quite  good.  It  would  appear  that  in  all  practical  cases  of  interest, 
the  approximate  pressure  solution  would  suffice. 

Figures  13  and  14  are  plots  of  the  pressure  fields  for  time  increments, 

At  =0.1  and  X  =  1  and  X  =  10.  Up  to  t  =  0.5,  the  pressure  distribution 

reflects  the  wave-propagation  nature  of  the  problem.  After  traversal  is 

completed,  i.e.,  t  >  0.5,  the  distribution  deviates  somewhat  from  the  usual 
experimental  observation  that  the  pressure  field  is  linear  for  any  given  t, 
e.g.,  Roberts  et  al.1® 

Figure  15  replots  the  pressure-time  variation  at  a  given  point  for 
X  =  0.1,  1  and  10.  The  approach  to  zero  pressure  at  the  center  of  the 
container  is  quite  abrupt  when  one  considers  the  amount  of  pressure  fluctuations 
still  present  at  x  =  0.4. 
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An  overview  of  Figures  11-15  would  indicate  that  the  Hayashi  approx¬ 
imation  for  the  pressure  is  excellent  for  X  >_  10.  To  correlate  the  pressure 
distribution,  however,  to  a  "maximum  acceleration",  which  is  due  entirely  to 
an  approximation,  is  not  self-consistent.  The  remarkable  agreement  with 
experimental  data  given  by  Hayashi,  reproduced  here  as  Figure  16,  should  be 
considered  fortuitious  at  this  time.  On  the  other  hand,  in  spite  of  its 
analytical  pathology,  the  concept  of  an  acceleration  beginning  at  zero  and 
in  phase  with  the  fluid  pressure  as  postulated  by  Hayashi  has  intuitive 
and  physical  appeal.  The  weaknesses  of  the  model  lie  obviously  in  neglecting 
the  container  mass  and  possibly  dissipation.  Typically,  one  can  add  a 
damper,  c,  in  parallel  with  the  spring,  k,  and  include  the  vessel  mass, 
i.e.,  the  boundary  condition  (18)  takes  the  dimensionless  form: 

u  ♦  nuT  ♦  xwuTT  *  *Ux(0,t)  -  ;x(1,t)}  ,  (40) 

where  u  is  the  container  to  fluid  mass  ratio  and  n  is  a  damping  parameter. 

In  a  preliminary  study,  the  author  has  shown  that  the  container  acceleration 
is  indeed  smoother  than  the  "jumps"  indicated  in  Figure  15,  e.g.,  the  acceler¬ 
ation  has  a  finite  rise  time  from  zero.  When  these  improvements  are  completed, 
the  model  first  proposed  by  Hayashi  might  indeed  still  make  the  discussion 
of  the  cavitation  model  very  much  easier  because  of  the  relatively  small 
number  of  dimensionless  parameters  present. 

TWO  DIMENSIONAL  MODELS 


The  use  of  a  fluid-filled  rigid  spherical  shell  under  impact  as 
a  model  of  craniocerebral  trauma  began  with  AnzeliuslS  and  GUttinger20. 
Phenomenologically,  this  model  is  not  very  different  from  the  one-dimen¬ 
sional  one  just  discussed.  Because  of  the  rigid  assumption,  a  blow 
at  one  point  (coup)  on  the  sphere  is  immediately  transmitted  to  the 
diametrically  opposite  pole  (contrecoup) .  Every  point  on  the  shell 
begins  to  transmit  energy  at  the  instant  of  impact.  Goldsmith^  sug¬ 
gested  the  obvious  necessity  of  modelling  the  skull  as  an  elastic 
shell.  Engin  and  Liu^1  justified  the  use  of  a  fluid-filled  elastic 
spherical  shell  as  a  model  for  head  injury  based  on  neuroanatomical 
and  analytical  considerations.  They  obtained  the  frequency  spectrum 
for  the  axisymmetric  free  vibration  of  the  fluid-solid  ensemble  and 
compared  the  results  for  the  cases  of  steel  and  water  and  skull  and 
water.  The  axisymmetric  response  of  a  thin,  elastic,  homogeneous  and 
completely  closed  spherical  shell  filled  with  an  inviscid  and  irrota- 
tional  fluid  subjected  to  a  radial  impulse  was  obtained  by  Engin22. 

The  point  of  view  taken  is  that  of  a  spherical  fluid  region  constrained 
by  very  complicated  boundary  conditions  consisting  of  the  two  shell 
differential  equations  (including  both  bending  and  membrane  effects) 
and  the  continuity  of  the  normal  velocity  at  the  fluid-solid  interface. 
Mathematically,  this  is  expressed  by  satisfying  the  wave  equation  in 
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The  pathological  nature  of  these  solutions  were  a  consequence  of 
the  mathematical  formulation  of  the  problem,  Wave  propagation  phenomena 
are  represented  by  solutions  of  hyperbolic  differential  equations  with 
prescribed  initial  and  boundary  data.  If  the  given  data  is  discontin¬ 
uous  (e.g.,  wave  motion  initiated  by  an  impulse:  uT(0)  >  V)  the  solu¬ 
tion  is  also  discontinuous.  The  function  f(x  t)  ♦  g(x  -  t)  is  a  gen¬ 
uine  solution  of  the  wave  equation  provided  f  and  g  are  twice  differen¬ 
tiable.  Otherwise,  the  function  ought  to  be  regarded  as  a  "solution 
in  the  general  sense"  or  a  weak  solution,  see  Courant  and  Hilbert. 18 

Because  of  the  inherent  weaknesses  of  the  above  one-dimensional 
model  due  to  both  the  mathematics  and  obvious  geometrical  restraints, 
further  improvement  of  this  model  appears  superfluous.  Typically,  if 
one  were  to  add  a  damping  element  in  parallel  or  in  series  with  the 
spring  and  include  the  inertia  of  the  container,  i.e.,  boundary  condi¬ 
tion  (18)  takes  the  form: 

u  +  nuT  ♦  viutt  ■  XUx(0,t)  -  Cx(1,t)}  ,  (40) 

where  n  and  u  are  the  dimensionless  damping  and  inertial  parameters 
respectively.  A  similar  analysis  as  given  earlier  would  yield  the 
presence  of  a  Dirac-delta  function  in  the  pressure  field,  which  is 
physically  untenable. 

It  is  indeed  regrettable  that  the  one-dimensional  continuum  model 
proposed  by  Hayashi  is  beset  by  these  mathematical  difficulties.  The 
model  has  the  advantage  of  simplicity.  If  it  had  satisfied  the  test 
of  self-consistency,  its  small  number  of  dimensionless  parameters  would 
have  made  discussions  of  the  cavitation  model  very  much  easier. 

TWO  DIMENSIONAL  MODELS 

The  use  of  a  fluid-filled  rigid  spherical  shell  under  impact  as  a 
model  of  craniocerebral  trauma  began  with  Anzeliusl®  and  Gtittinger20. 
Phenomenologically,  this  model  is  not  very  different  from  the  one-di¬ 
mensional  one  just  discussed.  Because  of  the  rigid  assumption,  a  blow 
at  one  point  (coup)  on  the  sphere  is  immediately  transmitted  to  the 
diametrically  opposite  pole  (contrecoup) .  Every  point  on  the  shell 
begins  to  transmit  energy  at  the  instant  of  impact.  Goldsmith*-*  sug¬ 
gested  the  obvious  necessity  of  modelling  the  skull  as  an  elastic 
shell.  Engin  and  Liu21  justified  the  use  of  a  fluid-filled  elastic 
spherical  shell  as  a  model  for  head  injury  based  on  neuroanatomical 
and  analytical  considerations.  They  obtained  the  frequency  spectrum 
for  the  axisymmetric  free  vibration  of  the  fluid-solid  ensemble  and 
compared  the  results  for  the  cases  of  steel  and  water  and  skull  and 
water.  The  axisymmetric  response  of  a  thin,  elastic,  homogeneous  and 
completely  closed  spherical  shell  filled  with  an  inviscid  and  irrota- 
tional  fluid  subjected  to  a  radial  impulse  was  obtained  by  Engin?2. 

The  point  of  view  taken  is  that  of  a  spherical  fluid  region  constrained 
by  very  complicated  boundary  conditions  consisting  of  the  two  shell 
differential  equations  (including  both  bending  and  membrane  effects) 
and  the  continuity  of  the  normal  velocity  at  the  fluid-solid  interface. 
Mathematically,  this  is  expressed  by  satisfying  the  wave  equation  in 
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axisymmetrical  spherical  coordinates  for  the  fluid  exactly,  i,e.,  the 
velocity  potential  of  the  fluid  in  the  Laplace-transformed  variable 
is  expanded  in  an  infinite  series  of  its  characteristic  functions. 

On  the  other  hand,  the  shell  equations  are  expanded  in  terms  of  Legendre 
and  associated  Legendre  polynomials,  which  are  the  characteristic  func¬ 
tions  for  the  in-vacuo  shell  differential  equations.  The  ensuing  in¬ 
finite  series  solutions  have  this  inherent  error  of  truncation. 

The  above  comments  have  a  familiar  ring  and  one  might  immediately 
and  legitimately  question  whether  or  not  the  solutions  thus  obtained 
fall  into  the  same  quagmire  as  the  infinite-series  solution  in  the 
Hayashi  model.  There  exists  no  exact  solution,  such  as  the  wave  propa¬ 
gation  result,  to  settle  the  question  definitively.  In  Engin's  solu¬ 
tion,  the  intracranial  pressure  required  the  evaluation  of  420  residues 
in  order  to  achieve  convergence  in  the  sense  of  Cauchy.  There  is  every 
physical  reason  to  believe  however,  that  these  solutions  will  not  con¬ 
tain  discontinuous  jumps.  By  comparison  with  the  one-dimensional  model, 
the  fluid-solid  model  has  the  following  advantages  in  this  regard: 

(a)  The  container  (shell)  is  elastic  and  hence  it  takes  a  finite 
time  for  the  stress  waves  to  travel  from  coup  to  contrecoup.  The  in- 
stantaneous  transmission  of  pressure  waves  due  to  the  assumption  of 
rigidity  of  the  container  is  one  contribution  to  the  jump  phenomena 
discussed  earlier  and  also  observed  in  the  solutions  due  to  GUttinger20. 

(b)  The  load  on  the  shell  is  spread  over  a  finite  area  (a  polar  cap) 
while  in  the  one  dimensional  model  it  is  confined  to  a  point,  a  physi¬ 
cally  untenable  state. 

(c)  The  obvious  geometrical  advantage  in  being  able  to  spread  the 
pressure  variation  as  a  function  of  the  polar  angle,  radius  and  time. 

The  results  of  Engin  confirmed  that  most  likely  locations  of  skull 
fracture  are  away  from  the  impact  polar  cap.  Specifically,  for  a  blow 
of  uniform  pressure  over  a  15°  cone  angle,  the  location  of  high  stresses 
on  the  skull  are  immediately  below  the  impact  area  and  in  a  ring  about 
45°  away  from  the  polar  axis.  The  second  interesting  result  is  the 
pressure  distribution  in  the  brain  during  a  closed-head  trauma.  If  one 
assumes,  for  the  moment,  that  injury  is  associated  with  negative  pres¬ 
sure,  then  the  possibility  of  intermediate  coup,  i.e.,  high  negative 
pressures  away  from  the  contrecoup  point,  has  been  predicted  by  the 
model.  Liu  and  Nelson2?  have  extended  Engin's  result  and  compared  the 
differences  in  pressure-wave  propagation  between  a  Dirac-delta  type 
impulse  to  a  more  realistic  finite  time  pulse.  Because  of  the  complexity 
of  the  wave  interaction  phenomena,  they  have  produced  a  "cartoon"  type 
movie  showing  the  progress  of  the  "isobars"  in  the  brain.  For  the  Dirac- 
delta  time  function,  the  frames  corresponding  to  contrecoup  and  inter¬ 
mediate  coup  are  shown  in  Figures  17  and  18.  In  these  figures,  the 
symbols  ♦,  0,  1,  2,  3,  4  and  5  represent  intervals  of  increasing  nega¬ 
tive  nondimensional  pressure  while  the  symbols  -,  »,  A,  B,  C,  D,  E,  and 
F  denote  increasing  positive  nondimensional  pressures.  The  dimension¬ 
less  pressure,  p^  and  time  t  are  related  to  the  physical  quantities  as 
follows: 

px  -  ps(l  -  v2)p/Epf  ,  (41) 
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t  ■  cst/a,  (42) 

where  ps  and  p^  are  the  densities  of  the  skull  and  fluid  respectively, 

E  and  v  are  the  Young's  modulus  and  Poisson's  ratio  of  the  shell  and 
p  is  the  magnitude  of  the  externally  applied  uniform  pressure,  cs  is 
the  wave  speed  in  the  shell  and  is  equal  to  {E/ps (1-v*) I5*,  a  is  the 
radius  of  the  shell  and  t,  the  time. 

We  note  from  Figures  17  and  13,  the  following: 

a)  There  exists  a  "ringing"  motion  in  the  fluid,  i.e.,  the  alter¬ 
nating  positive  and  negative  pressures  in  each  frame. 

b)  The  so-called  contrecoup  negative  pressure  is  low  in  magnitude 
(grade  1)  but  diffuse  over  a  fairly  large  area. 

c)  The  intermediate-coup  negative  pressure  is  high  in  magnitude 
(grade  S)  but  extremely  well-focused. 

The  response  of  the  fluid- sol id  interaction  system  to  a  more  realis¬ 
tic  finite  time-varying  pulse  is  obtained  by  convolution.  Using  the 
same  spatial  distribution  as  Engin,  but  instead  of  the  Dirac-delta 
function,  we  use 

T(t)  «  exp  (-4.73  x  10“6t2)  sin(w  x  10'3t)  (43) 

The  pressure  response  of  the  system  is  similar  except  that  the 
"ringing"  is  more  pronounced  and  the  contrecoup  area  less  diffused.  It 
was  evident  from  an  examination  of  the  frames  that  the  damage  criterion 
for  any  given  location  must  incorporate  the  following  factor:  the 
time  spent  beyond  a  critical  negative  pressure.  Physiological  consider¬ 
ations  collaborate  the  above  argument.  The  capillaries  will  not  burst 
even  for  a  large  transmural  pressure  if  the  negative  pressure  is  present 
for  only  a  very  short  time.  Finite  time  is  required  for  the  blood  to 
enlarge  the  capillary  walls  to  the  point  of  rupture.  On  the  other  hand, 
a  small  transmural  pressure  acting  over  a  relatively  long  period  of  time 
can  damage  the  vessel.  From  simple  strength  of  material,  the  static 
hoop  stress,  °,  for  a  thin-walled  tube  is 

a  *  Apd/2h,  (44) 

where  Ap  is  the  transmural  pressure,  h/d  is  the  wall  thickness  to  diame¬ 
ter  ratio.  According  to  Bv.rton24f  h/d  is  about  1/8  for  a  true  capillary. 
Its  thickness  is  practically  a  single  layer  of  endothelial  lining  cells. 
Very  little  force  is  required  to  deform  these  cells  and  they  are  believed 
to  play  a  very  minor  role  in  the  total  elasticity.  Without  going  into 
the  role  of  the  surrounding  tissue,  it  suffices  to  note  that  equation 
(44)  is  the  basis  of  those  who  claim  that  negative  pressure  is  the  damage 
criterion.  Liu  and  Nelson23  postulated  the  following  Cumraulative  Damage 
Criterion  for  closed-head  brain  injuries:  The  severity  of  the  brain 
injury  at  a  given  location  is  proportional  to  the  averaged  time  spent 
there  beyond  a  critical  negative  pressure,  p*.  This  criterion  is  il¬ 
lustrated  in  Figure  19.  The  negative  pressure  less  than  p*,  at  loca¬ 
tion  x0,  is  divided  into  equal  time  steps  At  apart.  Let  Pi(x0,At2), 
P2(xo»At2^*  Pi(xo»Ati)  be  tbe  lengths  of  a  series  of  chords 

representing  the  discrete  pressures  at  the  beginning  of  its  At^,  then 
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the  time-averaged  pressure  p(x0,t)  is  given  analytically  by: 

n  n 

P(x0,F)  *  £  p,(x  ,At.)H(p*  -  p,)/£At.,  (45) 

i»l  i»l 

where  H(p*  -  p^)  is  the  Heaviside  function,  i.e^_,  it_is  zero  if  p*  <  p^ 
and  is  unity  of  p*  >  p*.  The  maximum  value  of  p(x0,t)  is  assigned  a 
grade  of  5  and  a  decreasing  linear  scale  is  used  thereafter  for  the  lower 
pressures.  The  computations  were  terminated  at  the  end  of  the  applied 
pulse.  By  the  Cumulative  Damage  Criterion,  there  appears  two  foci  of  inter- 
mediate-coup  due  essentially  to  the  wave  propagation  and  interaction  phenomena 
and  con1  recoup  because  of  the  induced  accelerational  difference  between  the 
shell  a  d  the  fluid. 

The  situation  appears  paradoxical  at  this  point.  As  the  elasticity 
of  the  s.iell  is  taken  into  account,  the  analysis  seems  to  deviate  more 
and  more  from  the  comfortable  notion  of  a  linear  pressure  distribution 
reported  by  the  experimentalists,  see,  for  example,  the  works  of  Lingren^, 
Unterharnscheidt  and  SellierZb,  and  Roberts  et  al.16.  The  situation  is  not 
an  unfamiliar  one.  The  experimentalists  made  measurements  without  a  theoretical 
model,  or  at  best,  the  model  was  constructed  from  the  data  a  posteriori  without 
much  regard  for  the  fundamental  requirements  of  the  boundary-value  problem. 

On  the  other  hand,  the  analysts  have  not  been  concerned  with  the  experimental 
evidence  of  the  suitability  of  his  basic  assumptions  or  if  he  has,  found  the 
experimental  situation  too  difficult  to  formulate  in  analytical  terms.  Hope¬ 
fully,  this  disagreement  will  be  resolved  as  additional  experimental  and 
analytical  work  presently  in  progress  is  published. 

Recently,  Lee  and  Advani^7  partially  extended  the  model  of  Engin 
and  Liu^1  to  include  transverse  shear  and  rotational  inertia  in  the  shell 
approximation.  The  net  effect  of  this  improvement  is  to  include  wave¬ 
lengths  of  the  order  of  the  thickness  of  the  shell.  The  qualitative 
differences  between  their  result  and  Engin 's  were  negligible  for  radial 
displacement  and  fibre  stresses.  Unfortunately,  the  boundary-value  prob¬ 
lem  associated  with  the  fluid  pressure  determination  was  not  solved.  It 
is  precisely  in  the  fluid,  however,  where  the  "high  frequency"  effects 
would  be  most  manifest. 

Except  for  the  extremely  unusual  case  in  which  the  blow  is  delivered 
such  that  its  resultant  force  passes  through  the  center  of  mass  of  the 
head,  the  head  will  experience  both  a  translation  and  rotation.  In  the 
case  of  a  whiplash-type  of  injury,  one  also  expects  both  motions  to  occur. 

To  assess  the  relative  contributions  of  each  of  these  proposed  damage 
mechanisms  to  head  trauma,  Chan^8  is  presently  investigating  the  non- 
axisymmetric  response  of  a  fluid-filled  spherical  shell.  The  fluid  is 
considered  inviscid  and  irrotational  and  an  "exact"  six  mode  shell 
theory  is  used.  The  implications  of  the  latter  is  that  the  analysis 
will  be  valid  for  a  thick  shell,  i.e.,  its  radius  to  wall-thickness  is 
between  4  and  10.  The  loading  consists  of  an  axisymmetric  pulse  perpen¬ 
dicular  to  the  surface  and  a  traction  tangential  to  it. 
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ROTATIONAL  IMPACT  MODELS 


Because  of  the  very  small  shear  modulus  of  the  brain  compared  with 
its  bulk  modulus,  Holbourn29  was  the  first  to  postulate  the  possibility 
of  brain  damage  due  to  rotational  acceleration  alone.  Von  Gierke3  gave, 
for  soft  tissue,  a  value  for  volume  compressibility  of  2.6  x  1010  dyne/cm2 
as  compared  to  a  shear  elasticity  of  2.5  x  104  dyne/cm2.  He  further  noted 
that  when  soft  tisc  .e  is  subjected  to  a  blow  most  of  the  energy  spreads 
through  the  tissue  in  the  form  of  shear  waves  and  not  of  compression  waves. 
Only  very  rapid  blows  produced  compression  waves  of  appreciable  magni¬ 
tude  in  tissue.  The  implications  for  energy  transmission  in  the  soft 
brain  tissue  is  obvious.  Engin  and  Wang3^  gave  a  theoretical  analysis 
of  the  driving  point  impedance  test  for  a  continuum  model  of  the  brain. 
Relationships  were  derived  from  which  some  of  the  viscoelastic  para¬ 
meters  can  be  determined  from  the  experimental  data. 

31 

Holboum's  observation  has  been  pursued  further  by  Ommaya  et  al. 
and  Unterharnscheidt  and  Higgins3^.  The  mathematical  models  used  were 
the  single  degree  of  freedom  rotational  version  of  equation  (1)  with  the 
angular  head  displacement  0  as  the  dependent  variable  and  the  angular 
acceleration  as  input.  The  motion  of  the  head  is  modelled  as  a  rigid 
body  constrained  by  a  torsional  spring.  The  failure  criterion  is  a 
clinical  indicator  of  concussion.  The  resultant  data  is  correlated  with 
the  parameters  of  the  input  pulse  in  a  form  similar  to  Figure  3. 

Martinez  and  Garcia33  used  a  three -degrees -of -freedom  nonlinear 
model  to  describe  the  whiplash  phenomena.  The  head  is  represented  by 
a  mass  and  the  resistance  of  the  neck  to  rotation  and  shear  are  repre¬ 
sented  by  torsional  and  cantilevered  springs  respectively.  Roberts 
et  al.3*  divided  the  head  into  two  rigid-body  compartments:  the  skull 
and  the  brain.  A  torsional  spring  placed  at  the  center  of  the  head 
simulates  the  resistance  of  the  brain  to  rotation  relative  to  the  skull. 

The  resistance  to  linear  displacement  is  simulated  by  linear  spring 
connecting  the  brain  to  the  skull.  A  parametric  study  was  made  using 
many  assumed  input  system  constants.  Hayashi?5  noted  that  for  very 
short-duration  rotational  impact,  the  severest  shear  deformation  is  in¬ 
duced  in  the  boundary  layer  near  the  inner  surface  of  the  skull.  From 
this  observation,  he  constructed  a  model  consisting  of  two  concentric 
cylinders,  which  are  connected  elastically  in  shear  by  the  brain  materi¬ 
al.  The  cylinders  are  assumed  to  be  rigid  bodies:  the  outer  shell  rep¬ 
resenting  the  skull  and  the  inner  one  the  essentially  undeformed  portion 
of  the  brain.  Using  estimated  input  data  and  a  data  point  on  the  toler¬ 
ance  curve  due  to  Ommaya  et  al.31,  he  obtained  through  his  analysis  a 
tolerance  relationship: 

0*t2  -  2.52,  (46) 

where  0*  is  the  critical  or  maximum  input  angular  acceleration  and  tp  is 
the  impact  duration. 

The  first  two-dimensional  continuum  model  for  torsional  impact  to 
the  head  was  proposed  by  Lee  and  Advani36.  The  model  is  a  linear  elas¬ 
tic  sphere  subjected  to  a  step  angular  acceleration  about  a  diametrical 
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axis.  The  displacement  is  prescribed  to  be  zero  at  the  surface  of  the 
sphere.  The  results  are  expressions  for  the  dynamic  torsional  displace¬ 
ment  and  shear  stress.  The  maximum  torsional  displacement  occurs  ap¬ 
proximately  at  normalized  radius  r  ■  0.42.  There  is  a  rapid  buildup  of 
shear  stresses  at  the  surface  of  the  r-phere.  If  the  time  for  a  shear 
wave  to  travel  from  the  surface  to  the  center  is  normalized  to  unity, 
then  for  nondimensional  time  of  less  than  one  (t  <  1),  the  shear  stress 
decreases  with  radius,  i.e.,  reaching  zero  at  the  appropriate  point.  Of 
course  nothing  happens  prior  to  the  arrival  of  the  stress  wave.  For 
t  >  1,  the  stress  distribution  is  almost  linear  with  respect  to  the  zero 
at  the  origin.  In  addition,  the  linear  viscoelastic  response  of  the 
sphere  is  studied  by  using  a  superposition  principle.  The  numerical 
results  for  a  Maxwell  fluid  model  indicate  an  almost  logarithmic  atten¬ 
uation  of  shear  stress  with  respect  to  the  dimensionless  radius. 

OTHER  CRITICAL  ORGANS  AND  CONCLUDING  REMARKS 

While  attention  has  been  focused  on  the  situation  where  the  spine 
and  head  are  the  critical  organs,  all  the  perspectives  gained  hopefully 
have  implications  in  other  impact  situations.  For  instance,  the  criti¬ 
cal  organ  in  human  tolerance  to  blast  exposure  is  most  probably  the  lung. 
Sass37  has  contended  that  the  lung  is  the  most  susceptible  organ  to  in¬ 
jury  because  of  its  low  density  and  elasticity  and  non-rigid  envelope, 
while  the  brain  is  the  least  vulnerable  because  of  its  essential  incom¬ 
pressibility  and  rigid  enclosure.  Experimental  and  theoretical  support 
for  this  contention  has  not  been  conclusive.  Failure  criterion  in  ma¬ 
terials,  whether  biological  or  not,  is  given  in  terms  of  stresses,  which 
are  then  related  to  the  strains  through  constitutive  equations.  The 
strain-displacement  relationships  provide  the  final  link  to  the  relative 
motion  of  the  body  tissues.  Unless  and  until  the  intermediate  relation¬ 
ships  are  known,  excessive  displacement  or  rate  of  displacement  can 
not  by  itself  be  a  criterion  of  injury  susceptibility.  A  good  illus¬ 
trative  example  is  the  predominance  of  compression  fracture  of  the  ver¬ 
tebral  end-plates  aid  anterior  lip  fractures  in  caudocephalad  accelera¬ 
tion.  Without  doubt,  in  this  mode  of  acceleration,  relatively  large 
displacements  of  the  lungs  and  diaphragm  do  take  place  as  evidenced  by 
some  occurrances  of  lung  and  visceral  injury  in  pilot  ejection.  Of 
course,  these  critical  organs  can  be  analyzed  in  a  similar  way  as  the 
spine  and  head.  The  question  is  always  one  of  priority. 

A  synopsis  of  the  exposure  limits  for  human  tolerance  to  blast  will 
reveal  a  similar  hierarchy  of  idealizations  in  the  modelling  process. 

The  initial  linearized  model  due  to  von  Gierke-*  involved  the  concept  of 
equal  strain  in  lung  tissue.  Tolerance  curves  plotting  the  maximum  over¬ 
pressure  versus  the  duration  of  exposure  were  given.  An  improved  model 
due  to  Holladay  and  Bowen37  accounted  for  the  nonlinear  effects  due  to 
large  amplitude  pressure  waves  and  severe  compression  of  the  lung  volume. 
Further  refinement  of  the  model  probably  lies  in  the  continuum  notion 
that  there  exists  a  wave  transmission  through  the  skin,  the  neuromuscu¬ 
lature,  the  thoracic  cage  and  the  lung  tissue  while  simultaneously  the 
blast  wave  works  its  way  down  the  bronchial  tree.  The  possible  inter¬ 
action  of  these  two  waves  of  differing  speeds  might  be  reflected  in  the 
intrapulmonary  pressure  measurements. 
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Figure  2,  Typical  tolerance  graph  based  on  equal  tissue  strain  assump 
tion.  Aq  and  tj  are  the  peak  acceleration  and  duration  of  the  applied 
pulse  respectively.  tc  is  the  critical  pulse  duration. 
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Figure  3.  Mans  distribution  of  the  human  cadaver  trunk.  Dotted  line 
is  data  assumed  in  Ref.  9. 
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Figure  4.  Mass  eccentricity  distribution  of  the  human  cadaver  trunk. 
Dotted  line  is  data  assumed  in  Ref.  9. 
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Figure  5.  Mass  moments  of  inertia  distribution  of  the  human  cadaver 
trunk.  Dotted  line  is  data  assumed  in  Ref.  9. 
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Figure  6.  Moment  distribution  along  the  spine  due  to  a  lOg  caudocephalad 
(+GZ)  acceleration  pulse  at  60  and  90  milliseconds.  Dotted  lines 
denote  result  using  the  data  assumed  ii<  Ref.  9. 
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Figure  7.  Axial  force  distribution  along  the  spine  due  to  a  lOg  caudo 
cephalad  (+G2)  acceleration  pulse  for  different  times.  Dotted  lines 
denote  result  using  data  assumed  in  Ref.  9. 
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Figure  8.  Shear  force  distribution  along  the  spine  due  to  a  lOg 
caudocephalad  (+GZ)  acceleration  pulse  for  different  times.  Dotted 
lines  denote  result  using  data  assumed  in  Ref.  9. 
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Figure  9.  Time-history  of  configuration  changes  of  the  spine. 
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Figure  16.  Comparison  of  theory  with  experiment  according  to  Hayashi. 
PA*pressure  at  C0UPJ  Pg"pressure  at  contrecoup. 
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Figure  17.  Pressure  distribution  at  the  tine  of  approximate  contre 
coup.  Note  the  relatively  large  diffuse  area  at  grade  1. 


Figure  18.  Pressure  distribution  at  the  time  of  intermediate  coup 
Note  the  sharply*- focused  area  at  grades  5,  4,  3  and  2. 


Figure  19.  Cumulative  Damage  Criterion  for  a  closed  head  injury.  The 
damage  mechanism  is  proportional  to  the  averaged  time  spent  at  a  given 
location  beyond  a  critical  negative  pressure,  p*. 


